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We present the first results on the Dalitz-plot structure and improved measurements of the time- 
dependent CP-violation parameters of the process B° — > K%K%K% obtained using 468 x 10 6 BB 
decays collected with the BABAR detector at the PEP-II asymmetric-energy B factory at SLAC. The 
Dalitz-plot structure is probed by a time-integrated amplitude analysis that does not distinguish 
between B° and B° decays. We measure the total inclusive branching fraction B(B° —¥ KgKgKg) = 
(6.19 ± 0.48 ± 0.15 ± 0.12) x 10 -6 , where the first uncertainty is statistical, the second is systematic, 
and the third represents the Dalitz-plot signal model dependence. We also observe evidence for the 
intermediate resonant states /o(980), /o(1710), and /2(2010). Their respective product branching 
fractions are measured to be (2.70 t\H ± 0.36 ± 1.17) x 10~ 6 , (0.50 ± 0.04 ± 0.10) x 10~ 6 , 
and (0.54 ±£^±0.03 ±0.52) x 10~ 6 . Additionally, we determine the mixing-induced CP-violation 
parameters to be 5 = -0.94 j* ± 0.06 and C = -0.17 ± 0.18 ± 0.04, where the first uncertainty 
is statistical and the second is systematic. These values are in agreement with the standard model 
expectation. For the first time, we report evidence of CP violation in B° — > KgKgKg decays; CP 
conservation is excluded at 3.8 standard deviations including systematic uncertainties. 

PACS numbers: 13. 66. Be, 14.40.-n, 13.25. Gv, 13.25.Jx, 13.20.Jf 



I. INTRODUCTION 

Over the past ten years, the B factories have shown 
that the Cabibbo-Kobayashi-Maskawa paradigm in the 
standard model (SM), with a single weak phase in the 
quark mixing matrix, accounts for the observed CP- 
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symmetry violation in the quark sector. However, there 
may be other CP-violating sources beyond the SM. 
Charmless hadronic B decays, like B Q — > K^K^Kg, 
are of great interest because they are dominated by 
loop diagrams and are thus sensitive to new physics ef- 
fects at large energy scales [l|. In the SM, the mixing- 
induced CP-violation parameters in this decay are ex- 
pected to be the same, up to ~ 1% [2], as in the tree- 
diagram-dominated modes such as B° — > J/tpKg. Both 
BABAR [3j and Belle [4[ have previously performed time- 
dependent CP-violation measurements of the inclusive 
mode B° — > KgKgKg, which is permissible because the 
final state is CP definite 

The structure of the Dalitz plot (DP), however, is of 
interest; although the time-dependent CP- violation pa- 
rameters S and C [see Eq. ([33]) ] can be measured inclu- 
sively without taking into account the phase space, dif- 
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ferent resonant contributions may have different values 
of these parameters in the presence of new physics. The 
statistical precision is not sufficient to perform a time- 
dependent amplitude analysis, but as we show below, 
it is possible to extract branching fractions from reso- 
nant contributions to the decay using a time-integrated 
amplitude analysis. Additionally, the amplitude analy- 
sis could shed light on the controversial /x(1500) res- 
onance: recent measurements of B° — > K + K~ K% and 
B ± ->■ K+K-K± from BaBAR @-[| and Belle [| E3] 
have shown evidence of a wide structure in the m K + K - 
spectrum around 1.5 GeV. In these measurements, it was 
assumed that this structure is a single scalar resonance; 
however a vector hypothesis could not be ruled out. The 
BABAR measurement of B + — > K + K~n + [llj appears to 
show an enhancement around 1.5 GeV, while the BABAR 
analysis of B ± — > KgKgW^ [l~2| finds no evidence of a 
possible /x(1500), suggesting that the structure is cither 
a vector meson or something exotic. An amplitude anal- 
ysis of B° — > KgKgKg will provide further insight into 
the nature of this structure, as only intermediate states 
of even spin are permitted due to Bose-Einstein statis- 
tics; an observation of the /x(1500) decaying to K S K S 
would require an even-spin state. Finally, the amplitude 
analyses of B — > Ktttt and B — > KKK modes may be 
used to extract the Cabibbo-Kobayashi-Maskawa angle 

This paper presents the first amplitude analysis and 
the final BABAR update of the time-dependent CP- 
asymmetry measurement of B° — > KgKgKg using the 
full T(AS) dataset. The amplitude analysis is time- 
integrated CP-averaged (i.e., it does not use flavor- 
tagging information to distinguish between B° and B° 
mesons). It takes advantage of the interference pattern 
in the DP to measure relative magnitudes and phases for 
the different resonant modes using B° — ► K S K S K S de- 
cays with K% -> TT+TT-, denoted by B° -> 3K^(tt + tt-). 
The magnitudes and phases are then translated into indi- 
vidual branching fractions for the resonant modes. The 
time-dependent analysis extracts the S and C parame- 
ters by modeling the proper-time distribution. This part 
of the analysis uses both B° — > 3A"g(7r + 7r~) events and 
events where one of the Kg mesons decays to 7r°7r° , de- 
noted by B° -> 2Ag(7r+7r-)Ag(7r 7r°). 

In Sec. [II] we briefly describe the BABAR detector and 
the dataset. The amplitude analysis is described in 
Sec. |nT] and the time-dependent analysis in Sec. IIV1 Fi- 
nally we summarize the results in Sec. |Vl 



II. THE BABAR DETECTOR AND DATASET 

The data used in this analysis were collected with 
the BaBAR detector at the PEP-II asymmetric-energy 
e + e~ storage ring at SLAC. The sample consists of an 
integrated luminosity of 426.0 fb _1 , corresponding to 
(467.8 ± 5.1) x 10 6 BB pairs collected at the T(4S) res- 
onance ( "on-resonance" ) , and 44.5 fb _1 collected about 



40 MeV below the T(45) ("off-resonance"). 

A detailed description of the BABAR detector is pre- 
sented in Ref. [l4|. The tracking system used for track 
and vertex reconstruction has two components: a sil- 
icon vertex tracker and a drift chamber, both operat- 
ing within a 1.5 T magnetic field generated by a super- 
conducting solenoidal magnet. A detector of internally 
reflected Cherenkov light associates Cherenkov photons 
with tracks for particle identification. The energies of 
photons and electrons are determined from the measured 
light produced in electromagnetic showers inside a Csl 
crystal electromagnetic calorimeter. Muon candidates 
are identified with the use of the instrumented flux return 
of the solenoid. 



III. AMPLITUDE ANALYSIS 

In Sees. IIII Al and ITlI Bl we describe the DP formalism 
and introduce the signal parameters that are extracted 
from data. In Sec. IIII CI we describe the requirements 
used to select the signal candidates and suppress back- 
grounds. In Sec. IIII Dl we describe the fit method and the 
approach used to account for experimental effects such as 
resolution. In Sec. IIII El we present the results of the fit, 
and finally, in Sec. IIII Fl we discuss systematic uncertain- 
ties in the results. 



A. Decay amplitudes 

The B° — > KgKgKg decay contains three identical 
particles in the final state and therefore the amplitude 
needs to be symmetrized. We consider the decay of a 
spin-zero B° into three daughters, Ka{l), Kg(2), and 
Kg(3), with four-momenta pi, p 2 , and p 3 . The decay 
amplitude is given by Q 

A[B° -+ A°(1)A°(2)A°(3)] (1) 

= {AilB ^ A°(1)A°(2)A°(3)] 

+ A 2 [B"^ A°(2)A°(3)A°(1)] 
+ A 3 [B°^ A°(3)A°(1)A°(2)]}, 

which takes into account the three permitted paths from 
the initial state to the final state. For instance for the 
B° decay this consists of an intermediate state K°K a K°. 
Since the labeling of the three identical particles is arbi- 
trary, we classify the final-state particles according to the 
square of the invariant mass, Sij, defined as 

Sij = Sji = ™>K%(i)K°(j) = (Pi+Pj) 2 ' ( 2 ) 

where i and j are the Kg indices. We use as independent 
(Mandelstam) variables the minimum and the maximum 



7 



of the squared masses s n 



= min(s 12 ,S23,Si3), 
= max(si2,s 2 3,si 3 )- 



(3) 



The third (median) invariant squared mass s mo d can be 
obtained from energy and momentum conservation: 



(4) 



The differential B meson decay width with respect to the 
variables defined in Eq. ([3]) (i.e., the DP variables) reads 



dT(B -+ K S K S K° S ) 



1 



(2tt) 3 32m 3 BO 



^min^max ) 



(5) 



where A is the Lorentz-invariant amplitude of the three- 
body decay. This amplitude analysis does not take into 
account any flavor tagging or time dependence, thus it 
is CP averaged and time integrated. The term \A\ 2 is 
therefore simply the average of squares of the contribu- 
tions A[B° K S K S K S ] and A[B° K%K%K%]. 

The choice of the variables s m ; n and s max gives a 
uniquely defined coordinate in the symmetrized DP. 
Therefore only one sixth of the DP is populated, i.e., 
the event density is six times larger compared to an am- 
plitude analysis involving three distinct particles. 

We describe the distribution of signal events in the 
DP using an isobar approximation, which models the to- 
tal amplitude as resulting from a coherent sum of am- 
plitudes from the N individual decay channels of the B 
meson, either into an intermediate resonance and a bach- 
elor particle or in a nonresonant manner: 



A(s 



mm ; J max J 



N 

E 



min ; '-'max ) 



(6) 



Here Fj (described in detail below) are DP-dependent 
amplitudes containing the decay dynamics and Cj are 
complex coefficients describing the relative magnitudes 
and phases of the different decay channels. This descrip- 
tion, which contains a single complex number Cj for each 
decay channel regardless of the B flavor (B° or B°), re- 
flects the assumptions of no direct CP violation and of 
a common weak phase for all the decay channels. With 
this description we cannot extract any weak phase in- 
formation; this would require using per-i? flavor complex 
amplitudes. The factor Fj contains strong dynamics only, 
and thus does not change under CP conjugation. 

Intermediate resonances decay to K°K°. In terms of 
the isobar approximation, the amplitude in Eq. (fTJ) for a 
resonant state j becomes 



A[B° -> K° s (l)K%(2)K° 8 (3)] 



(7) 



OC Cj [Fj(su,Si 3 ) + Fj(si2,S 2 3) +-P)(S13,S23)] • 

This reflects the fact that it is impossible to associate a 
given Kg to a flavor eigenstate K° or K°. In practice, 



this sum of three Fj terms, corresponding to an even- 
spin resonance, is implicitly taken into account by the 
description in terms of s m i n and s max . 

The Fj terms are represented by the product of the 
invariant mass and angular distributions; i.e., 

^■(smin,Sma X ,i) = Rj (m)X L ( \p * \ r')X L ( | q I r)T 3 (L, p, q ) , 

(8) 

where 

(i) m is the invariant mass of the decay products of the 
resonance; 

(ii) Rj (m) is the resonance mass term or "line shape" , 
e.g., relativistic Breit-Wigner (RBW); 

(iii) L is the orbital angular momentum between the res- 
onance and the bachelor particle; 

(iv) p * is the momentum of the bachelor particle evalu- 
ated in the rest frame of the B; 

(v) p and q are the momenta of the bachelor particle 
and one of the resonance daughters, respectively, 
both evaluated in the rest frame of the resonance; 

(vi) X L (\p*\ r') and X L (\q\ r) are Blatt-Weisskopf bar- 
rier factors [lH with barrier radii of r and r' , and 



(vii) Tj(L,p,q) is the angular distribution: 
L = : Tj = I, 
L = 2 : T, 



(9) 



,-3 [3(p-qf-(\p\\q\) 2 }- (10) 



The Blatt-Weisskopf barrier factor is unity for all the 
zero-spin resonances. In our analysis it is relevant only 
for the /2(20I0). Since for this resonance r and r' are not 
measured, we take them both to be 1.5 GeV -1 and vary 
by ±0.5 GeV" 1 to estimate the systematic uncertainty. 

The helicity angle of a resonance is defined as the an- 
gle between p and q. Explicitly, the helicity angle 9 
for a given resonance is defined between the momenta 
of the bachelor particle and one of the daughters of the 
resonance in the resonance rest frame. Because of the 
identical final-state particles this definition is ambiguous, 
but the ambiguity disappears because of the description 
of the DP in terms of s m - m and s max . There are three 
possible invariant-mass combinations: s m ; n , s mo d, and 
Smax- We denote the corresponding helicity angles as 
0minj $med, and max - The three angles are defined be- 
tween and 7r/2. 

As the present study is the first amplitude analysis of 
this decay, we use the method outlined in Sec. MI D 31 
to determine the contributing intermediate states. The 
components of the nominal signal model are summarized 
in Table |U 

For most resonances in this analysis the Rj are taken 
to be RBW [16| line shapes: 



Rj(m) = 



(m§ — m 2 ) — imoT(m) ' 



(11) 
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where mo is the nominal mass of the resonance and T(m) 
is the mass-dependent width. In the general case of a 
spin- J resonance, the latter can be expressed as 



T(m) 



2J+1 



A 



fm \ 

Vm J XjQqhlry 



(12) 



The symbol To denotes the nominal width of the reso- 
nance. The values of too an d To are listed in Table HI 
The symbol go denotes the value of q when m = mo- 

For the /o(980) line shape the Flatte form [13] is used. 
In this case the mass-dependent width is given by the 
sum of the widths in the 7T7T and KK systems: 



where 



T(m) = r w7r (m) + T KK (m), 



r 7r7r (m) = g n ( i sjl - im^o/m 2 + 



(13) 



(14) 



-^/l-4m2 ± /m 2 , 



9k I - -im 2 K± /m 2 + 



(15) 



-^l-Am 2 K0 /m 2 



The fractional coefficients arise from isospin conservation 
and g n and gx are coupling constants for which the values 
are given in Table HI The nonresonant (NR) component 
is modeled using an exponential function: 



therefore define /i m in and h max as cos 6> m ; n and cos 9 n 
respectively, and apply the transformation 



(s 



mini i 1 1 : ' s i 



ij : in? ^max 



(17) 



The (ft-min, /imax) plane is referred to as the square Dalitz 
plot (SDP), where both /i m ; n and /i max range between 
and 1 due to the convention adopted for the helicity 
angles (see Fig.rjJ. Explicitly, the transformation is 



mm V'-- max 



^mcd ) 



•Smin^ 4t?1^o ^min 



(18) 



1 



(^^o ^min)^ 
Smax(^mcd ^min) 

— X 

l K%' 



l K°- 



(19) 



^(m 2 BO -m 2 K o 



4to 



2 i 



where the numerators may easily be expressed in terms 
of s m ; n and s max using Eq. (QJ. The differential surface 
elements of the DP and the SDP are related by 



| det J | dh min d/i max , 



(20) 



where J = J (h mm , h mllx ) is the appropriate Jacobian 
matrix. The backward transformations s m i n (/i m i n , /imax) 
and s max (/i m in, /imax), and therefore the Jacobian 
| det J |, cannot be found analytically; they are obtained 
numerically. The variables /i m in and h max as a function 
of the invariant masses are shown in Fig. [T] together with 
the Jacobian. 



-Rnr(w) = e c 



(16) 



As in the resonant case, here m is the invariant mass 
of the relevant K^Kg pair. The parameter a is taken 
from the BABAR B + K + K~K + analysis @, H| and is 
given in Table HI This value was found to be compatible 
with the one resulting from varying a in the maximum- 
likelihood fit in the present analysis. There is no satisfac- 
tory theoretical description of the NR component; it has 
to be determined empirically. The exponential function 
of Eq. (TlT)f was used by other amplitude analyses of B- 
meson decays to three kaons [a-llOj. Adopting the same 
parametrization for the NR term allows the comparison 
of results for other components. 



B. The square Dalitz plot 

We use two-dimensional histograms to describe the 
phase-space dependent reconstruction efficiency and to 
model the background over the DP. When the phase- 
space boundaries of the DP do not coincide with the his- 
togram bin boundaries this may introduce biases. We 



C. Event selection and backgrounds 



We reconstruct B° 



KgKgKg candidates from three 



Kg — > 7r + 7r~ candidates that form a good quality vertex; 
i.e., the fit of the B° vertex is required to converge and 
the x 2 probability of each K s vertex fit has to be greater 
than 10 -6 . Each K s candidate must have 7r + 7r~ invari- 
ant mass within 12.1 MeV/c 2 of the nominal K° mass [1 61 ] . 
and decay length with respect to the B vertex between 
0.22 and 45 cm. The last criterion ensures that the de- 
cay vertices of the B° and the K s are well separated. 
In addition, combinatorial background is suppressed by 
selecting events for which the angle between the momen- 
tum vector of each Kg candidate and the vector con- 
necting the beamspot and the K s vertex is smaller than 
0.0185 radians. We ensure a good B vertex fit quality by 
requiring that the charged pions of at least one of the Kg 
candidates have hits in the two inner layers of the vertex 
tracker. 

A B meson candidate is characterized kine- 
matically by the energy-substituted mass toes = 
\J (s/2 + pi ■ pb) 2 /E 2 — p 2 B and the energy difference 
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TABLE I: Parameters of the DP model used in the fit. The Blatt-Weisskopf barrier parameters (r and r') of the /2(2010), 
which have not been measured, are varied by ±0.5 GeV - for the model uncertainty. 



Resonance 


Parameters 


Line shape 


Reference 


/o(980) 


m = (965 ± 10) MeV/c 2 
= (165 ± 18) MeV/c 2 
g K = (695 ± 93) MeV/c 2 


Flatte 
Eq. dHJ) 


nsi 


/o(mo) 


m = (1724 ±7) MeV/c 2 
To = (137 ± 8) MeV/c 2 


RBW 
Eq. fllj 


[16] 


h (2010) 


mo = (201ll^)MeV/c 2 
r = (202 ± 60) MeV/c 2 
r = r' = 1.5 GeV" 1 


RBW 
Eq. QJ 


[16] 


NR decays 


a = (-0.14 ± 0.02) GeV~ 2 c 4 


Exponential NR 
Eq. (IB) 


[8J 


XcO 


m = (3414.75 ± 0.31) MeV/c 2 
To = (10.2 ± 0.7) MeV/c 2 


RBW 
Eq. (JTTJ) 


rie] 




s max [GeV 2 /c 4 ] h min 

FIG. 1: (color online). Lines of constant helicity angle in the Dalitz plot of s m j n versus s max (left), and the magnitude of the 
Jacobian (gray scale on the right) mapping (s m i n , s max ) to (/i m i n , ftmax). For the latter see Eq. (|20[) . 



AE = Eg — \\fs, where (E b ,Pb) and (Ei,pi) are the 
four- vectors in the laboratory frame of the i?-candidate 
and the initial electron-positron system, respectively, 
and pb is the magnitude of pu- The asterisk denotes 
the T(45) frame, and s is the square of the invariant 
mass of the electron-positron system. We require 
5.27 < m ES < 5.29GeV/c 2 and \ AE\ < 0.1 GeV. Follow- 
ing the calculation of these kinematic variables, each of 
the B candidates is refitted with its mass constrained 
to the world average value of the B meson mass [IB] in 
order to improve the DP position resolution, and ensure 
that Eq. (|4]) holds. The sideband used for background 
studies is in the range 5.20 < to E s < 5.27GeV/c 2 and 
\AE\ < 0.1 GeV. 

Backgrounds arise primarily from random combina- 
tions in continuum e + e~ — ¥ qq events (q = d,u, s, c). 
To enhance discrimination between signal and contin- 



uum background, we use a neural network (NN) [19| to 
combine four discriminating variables: the angles with 
respect to the beam axis of the B momentum and B 
thrust axis in the T(AS) frame, and the zeroth- and 
second-order monomials Lo,2 of the energy flow about 
the B thrust axis. The monomials are defined by L n = 
J2iPi x |cos#j| n , where 8i is the angle with respect to 
the B thrust axis of track or neutral cluster i and pi 
is the magnitude of its momentum. The sum excludes 
the B candidate and all quantities are calculated in the 
T(AS) frame. The NN is trained with off-resonance data, 
sideband data and simulated signal events that pass the 
selection criteria. Approximately 0.5% of events passing 
the full selection have more than one candidate. When 
this occurs, we select the candidate for which the error- 
weighted average of the masses of the K% candidates is 
closest to the world average K% mass hj. With the 
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above selection criteria, we obtain a signal reconstruc- 
tion efficiency of 6.6% that has been determined from 
a signal Monte Carlo (MC) sample generated using the 
same DP model and parameters as obtained from the 
data fit results. We estimate from this MC that 1.4% of 
the selected signal events are misreconstructed, and as- 
sign a systematic uncertainty (see Sec. lIIIFl) . We use MC 
events to study the background from other B decays (B 
background). We expect fewer than 6 such events in our 
data sample. As these events are wrongly reconstructed, 
the toes and AE distributions are continuumlike and as 
a result the events are mostly absorbed in the continuum 
background category. We assign a systematic uncertainty 
for B background contamination in the signal. 



D. The maximum-likelihood fit 

We perform an unbinned extended maximum- 
likelihood fit to extract the B° -> K S K° S K S event yield, 
as well as the resonant and nonresonant amplitudes. The 
fit for the amplitude analysis uses the variables toes and 
AE, the NN output, and the SDP variables to discrimi- 
nate signal from background. The selected on-resonance 
data sample is assumed to consist of signal and con- 
tinuum background. The feed-through from B decays 
other than the signal is found to be negligible. Misrecon- 
structed signal events are not considered as a separate 
event species, but are taken into account as a part of the 
signal. The likelihood function Ci for event i is the sum 



d = 



3 



Autoes, A£,NN, 

), (21) 



where j stands for the species (signal, continuum back- 
ground) and Nj is the corresponding yield. Each prob- 
ability density function (PDF) is the product of four 
individual PDFs: 

V) = V;(m ES ) V}{AE) P](NN) P}{h min ,h max ). (22) 

A study with fully reconstructed MC samples shows that 
correlations between the PDF variables are small and 
therefore we neglect them. However, possible small dis- 
crepancies in the fit results due to these correlations 
are accounted for in the systematic uncertainty (see 
Sec. lIIlFf . The total likelihood is given by 



The mo parameters for both toes and AE are free in 
the fit to data, while the other parameters are fixed to 
values determined from a fit to MC simulation. For the 
NN distributions of signal we use a histogram PDF from 
MC simulation. 

For continuum events the toes and AE PDFs are 
parametrized by an ARGUS shape function [20( and a 
straight line, respectively. The NN PDF is described by 
a sum of power functions: 

E(x; ci, a, b , h, b 2 , h, c 2 , e 3 ) = (25) 
cos 2 (ci) [cos 2 (a) Af(b , h) x b ° (1 - x) bl 
+ sm 2 (a)M(b 2 ,b 3 )x b ^ (l-x) 63 ] 
+ sin 2 (ci)AA(c 2 , C3 )x C2 (l-x) C3 , 

where x = (NN - NN min ) / (NN max - NN min ) and the J\f 
are normalization factors, computed analytically using 
the standard T function: 



Af(a,(3) = 



r(/3 + 2 + q) 
T(a + l)T(j3 + 1) ' 



(26) 



The parameters for all the continuum PDFs are deter- 
mined by a fit to sideband data and then fixed for the fit 
in the signal region. 



2. Dalitz-plot PDFs 

The SDP PDF for continuum background is a his- 
togram obtained from jjjes sideband on-resonance events. 
The SDP signal PDFs require as input the DP-dependent 
selection efficiency, e — £(h m i n ,h maxi ), that is described 
by a histogram and is taken from MC simulation. For 
each event we define the SDP signal PDF: 



^sig (^min ; ^max) (X £\hxnini ^max 



|-^(^min7 ^max 



C = exp(- 



E 



(23) 



The normalization of the PDF is implemented by nu- 
merical integration. To describe the experimental reso- 
lution in the SDP variables, we use an ensemble of two- 
dimensional histograms that represents the probability to 
reconstruct at the coordinate (h m - m \ /i ma x') an event that 
has the true coordinate (/i m i n , h max ). These histograms 
are taken from MC simulation and are convolved with 
the signal PDF. 



1. The m-ES, AE, and NN PDFs 

The toes and AE distributions of signal events are 
parametrized by an asymmetric Gaussian with power- 
law tails: 



Cr(x; m , a t , oy, a>i,a r ) 
(x - m ) 2 



(24) 



exp 



x — too < : i = I 
2of + aa(x - m ) 2 y\a;-TO O >0: i = r 



3. Determination of the signal Dalitz-plot model 

Using on-resonance data, we determine a nominal sig- 
nal DP model by making likelihood scans with various 
combinations of isobars. We start from a baseline model 
that includes / (980), XcO, and NR components. We 
then add another scalar resonance described by the RBW 
parametrization. We scan the likelihood by fixing the 
width and mass of this additional resonance at several 
consecutive values, for each of which the fit to the data 
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is repeated. All isobar magnitudes and phases are float- 
ing in these fits. From the scans we observe a significant 
improvement of the fit around a width and mass that are 
compatible with the values of the /o(1710) resonance 
After adding the /o(1710) to the nominal model we re- 
peat the same procedure for an additional tensor particle. 
We find that the /2(2010) has a significant contribution. 
The results of the likelihood scans are shown in Fig. [5] 
in terms of — 2Aln£ = — 21n£ — (— 21n£) m i n , where 
(— 21n£) m j n corresponds to the minimal value obtained 
in the particular scan. To conclude the search for pos- 
sible resonant contributions we add all well established 
resonances [l6| and check if the likelihood increases. We 
do not find any other significant resonant contribution 
but as we cannot exclude small contributions from the 
/o(1370), / 2 (1270), /a(1525), ao(1450), and / (1500) res- 
onances, we assign model uncertainties (see Sec. MI F[) 
due to not taking these resonances into account. 



E. Results 

The maximum-likelihood fit of 505 candidates results 
in a B° -> K%K%K% event yield of 200± 15 and a contin- 
uum yield of 305 ± 18, where the uncertainties are statis- 
tical only. The symmetrized and square Dalitz plots of a 
signal DP-model MC sample generated with the result of 
the fit to data are shown in Fig.[3J Figure 0] shows plots of 
AE, toes, and the NN for isolated signal and continuum 
background events obtained by the s Vlots [2l[ technique. 
Figure [5] shows projections of the data onto the invariant 
masses s min and s max . 

When the fit is repeated with initial parameter values 
randomly chosen within wide ranges above and below the 
nominal values for the magnitudes and within the [—it, 7t] 
interval for the phases, we observe convergence towards 
two solutions with minimum values of the negative log 
likelihood function — 21n£ separated by 3.25 units. In 
the following, we refer to them as Solution 1 (the global 
minimum) and Solution 2 (a local minimum). No other 
local minima were found. 

In the fit, we measure directly the relative magni- 
tudes and phases of the different components of the signal 
model. The magnitude and phase of the NR amplitude 
are fixed to 1 and 0, respectively, as a reference. In Fig. [5] 
we show likelihood scans of the isobar magnitudes and 
phases of all the resonances, where both solutions can be 
noticed. Each of these scans is obtained by fixing the cor- 
responding isobar parameter at several consecutive val- 
ues, for each of which the fit to the data is repeated. The 
measured relative amplitudes are used to extract the 
fit fraction defined as 

t'W = 7TB F *\ ' (.2°) 

c IJ. c v\ 1 V 1 ul 

where k, which varies from 1 to 5, represents an interme- 
diate state. Each fit fraction is a sum of three identical 



contributions, one for each pair of K s . The indices \x and 
v run from 1 to 15, as each of the five resonances con- 
tributes to three pairs of Kg, which correspond to the 
three terms (3fc — 2, 3fc — 1, and 3k) in each sum in the 
numerator of Eq. ([28]) . The dynamical amplitudes F are 
defined in Sec. MI Al and the terms 



F F*.7s • 

J - fi A y u " , min u " , max 



(29) 



are obtained by integration over the DP. The total fit 
fraction is defined as the algebraic sum of all fit fractions. 
This quantity is not necessarily unity due to the potential 
presence of net constructive or destructive interference. 

In order to estimate the statistical significance of each 
resonance, we evaluate the difference Aln£ between the 
log-likelihood of the nominal fit and that of a fit where 
the magnitude of the amplitude of the resonance is set 
to (this difference can be directly read from the like- 
lihood scans as a function of magnitudes in Fig. [5]) . In 
this case the phase of the resonance becomes meaning- 
less, and we therefore account for two degrees of freedom 
removed from the fit. The value 2 A ln£ is used to evalu- 
ate the p- value for 2 degrees of freedom; we determine the 
equivalent one-dimensional significance from this p- value. 

The results for the phase and the fit fraction are given 
in Table [TT1 for the two solutions; the change in likelihood 
when the amplitude of the resonance is set to and the 
resulting statistical significance of each resonance is given 
for Solution 1. 

As the fit fractions are not parameters of the PDF it- 
self, their statistical errors are obtained from the 68.3% 
coverage intervals of the fit-fraction distributions ob- 
tained from a large number of pseudoexperiments gen- 
erated with the corresponding solution (1 or 2). As ob- 
served in other three-kaon modes 043' the total FF 
significantly exceeds unity. 

In Table |TT] it can be seen that the two solutions differ 
mostly in the fraction assigned to the NR and the /o (980) 
components. Solution 1 corresponds to a small FF of the 
/o(980) and a large value for the NR, and Solution 2 has 
a large /o(980) fraction and a smaller NR fraction. Other 
three-kaon modes [6-10] favor the behavior of Solution 1. 



Generalizing Eq. (|28|) . we obtain the interference frac- 
tions among the intermediate decay modes k and j 



FF(k,j) 



3fc 

iU=3fe-2 



3j 



£2 



=3j-2 



(30) 



which are given in Table UTT1 for Solution 1. Unlike the 
total FF defined above, the elements of this matrix sum 
to unity. The large destructive interference between the 
io(980)Alj and the NR components appears clearly in the 
table. This is possible due to the large overlap in phase 
space between the exponential NR term and the broad 
tail of the /o(980) resonance above the KK threshold. 

Using the relative fit fractions, we calculate the branch- 
ing fraction B for the intermediate mode k as 



FF(k) x B(B° 



K%K%K%) 



(31) 
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FIG. 2: Two-dimensional scans of — 2Aln£ (gray scale) as a function of the mass and the width of an additional resonance. 
These scans were performed to look for an additional scalar resonance (left) and an additional tensor resonance (right). The 
baseline model of the scans for additional scalar resonances contains /o(980), XcO, and NR intermediate states. The baseline 
model of the scans for additional tensor resonances contains /o(980), XcO, NR, and /o(1710) intermediate states. The ellipses 
indicate the world average parameters [16|] for the /o(1710) and /2(2010) resonances that are added to the model. 
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FIG. 3: Symmetrized (left) and square (right) DP for MC simulated signal events using the amplitudes obtained from the fit 
to data. The low population in bins along the edge of the symmetrized DP is due to the fact that the phase space boundaries 
do not coincide with the histogram bin boundaries. 



where B(B° 
ing fraction: 



Kg KgKg) is the total inclusive branch- 



B(B° -+ K S K S K° S ) 



A si 



£N BS 



(32) 



We estimate the average efficiency e = 6.6% using a fully 
reconstructed DP-model MC sample generated with the 
parameters found in data. The results of the branching 
fraction measurements are shown in Table llVl As a cross 
check we attempt to compare our measured branching 
fractions to results from other measurements; however, 
many of the branching fractions for the decay into kaons 
of the resonances included in our model are not (or are 
only poorly) measured (marked as "seen" in Ref. 16]). 



An exception is the charmonium state x c o, for which the 
measured value is B(xco -> K S K%) = (3.16 ± 0.18) x 
10 -3 [l6j]. We can then use the BaBAR measurement of 
B(B° XcoK ) = (142 ±g ± 8 ± 16 ± 12) x 10" 6 H3 
to calculate B[B° -> X co(-^ K S K^)K S ] = ±B(B° -> 
XcoA' ) x B{ X co -> A° K° s ) = (0.224 ± 0.078) x lO" 6 , 
which is consistent with our measured branching fraction, 
given in Table IIV1 

An interesting conclusion from this first amplitude 
analysis of the B° — > KgKgKg decay mode is that we do 
not need to include a broad scalar fy (1500) resonance, as 
has been done in other measurements |6T-fl0l|. to describe 
the data. The peak in the invariant mass between 1.5 
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FIG. 4: (color online). a Vlots (points with error bars) and PDFs (histograms) of the discriminating variables: toes (top), AE 
(middle), and NN (bottom), for signal events (left) and continuum events (right). Below each bin are shown the residuals, 
normalized in error units. The horizontal dotted and full lines mark the one- and two-standard deviation levels, respectively. 



and 1.6 GeV/ c 2 can be described by the interference be- 
tween the /o(1710) resonance and the nonresonant com- 
ponent. However minor contributions from the / 2 (1525) 
and /o(1500) resonances to this structure cannot be ex- 
cluded. 



F. Systematic uncertainties 

Systematic effects are divided into model and exper- 
imental uncertainties. Details on how they have been 
estimated are given below and the associated numerical 
values are summarized in Table fVl 



1. Model uncertainties 

We vary the mass, width, and any other parameter of 
all isobar fit components within their errors, as quoted 
in Table HI and assign the observed differences in our 
observables as the first part of the model uncertainty 
("Model" in Table E]). To estimate the contribution to 
B° — > KgKgKg from resonances that are not included 
in our signal model but cannot be excluded statistically, 
namely the / (1370), / 2 (1270), / 2 (1525), a (1450), and 
/o(1500) resonances, we perform fits to pseudoexperi- 
ments that include these resonances. The masses and 
the widths are taken from [l6|, except for the /o(1370) 
for which we take the values from [23J . We generate pseu- 
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FIG. 5: (color online). Projections onto ^/s m i n (left) and ^/s m „ (right). On-resonance data are shown as points with error 
bars while the dashed (dotted) histogram represents the signal (continuum) component. The solid-line histogram is the total 
PDF. Below each bin are shown the residuals, normalized in error units. The horizontal dotted and full lines mark the one- 
and two-standard deviation levels, respectively. 



TABLE II: Summary of measurements of the quasi-two-body 
parameters. The quoted uncertainties are statistical only. 
The change in the log-likelihood (— 2Aln£) corresponds to 
the case where the magnitude of the amplitude of the reso- 
nance is set to 0. This number is used for the estimation of 
the statistical significance of each resonance. 

Mode Parameter Solution 1 Solution 2 

/o(980)A° FF 0.44+°?° 1.03 1°;^ 

Phase [rad] 0.09 ±0.16 1.26 ±0.17 

-2Aln£ 11.7 

Significance [a] 3.0 

/o(1710)Ag FF 0.07±°:£ 0.09 1°;°* 

Phase [rad] 1.11 ±0.23 0.36 ± 0.20 

-2Aln£ 14.2 

Significance [a] 3.3 

f 2(2010) Ks FF 0.09t°;o3 0.10 ±0.02 



Phase [rad] 2.50 ± 0.20 1.58 ± 0.22 
-2Aln£ 14.0 
Significance [a] 3.3 



NR 


FF 


2 lfi+ - 36 
z - lu -0.37 


1 oy+0.26 
-0.21 




Phase [rad] 


0.0 


0.0 




-2Aln£ 


68.1 






Significance [a] 


8.0 




XcoAg 


FF 


07+ 04 
u - u ' -0.02 


0.07 ±0.02 




Phase [rad] 


0.63 ± 0.47 


-0.24 ±0.52 




-2Aln£ 


18.5 






Significance [a] 


3.9 






Total FF 


2-84 1°;£ 


2.66 t°rjf 



doexperiments with the additional resonances, where the 
isobar magnitudes and phases have been determined in 
fits to data, and fit these datasets with the nominal 
model. We assign the induced shift in the observables 
as a second part of the model uncertainty. 

2. Experimental systematic uncertainties 

To validate the analysis procedure, we perform fits on 
a large number of pseudoexperiments generated with the 
measured yields of signal events and continuum back- 
ground. The signal events are taken from fully recon- 
structed MC that has been generated with the fit result 
to data. We observe small biases in the isobar magni- 
tudes and phases. We correct for these biases by shift- 
ing the values of the parameters and assign to this pro- 
cedure a systematic uncertainty, which corresponds to 
half the correction combined in quadrature with its er- 
ror. This uncertainty accounts also for correlations be- 
tween the signal variables, wrongly reconstructed events, 
and effects due to the limited sample size ( "Fit Bias" in 
Table E| . 

From MC we estimate that there are six B background 
events in our data sample. To determine the bias in- 
troduced by these events, we add B background events 
from MC to our data sample, and fit it with the nominal 
model. We then assign the observed differences in the 
observables as a systematic uncertainty ( "i?-bkg" in Ta- 
ble EJ . We assign a systematic uncertainty for all fixed 
PDF parameters by varying them within their uncertain- 
ties according to the covariance matrix. 

We vary the histogram PDFs, i.e., the SDP PDF for 
continuum and the NN PDF for signal ("Discr. Vars" 
in Table El). The to ES dependence of the SDP PDF for 
continuum was found to be negligible. We account for 
differences between simulation and data observed in the 
control sample B° -t J/ipK° s ("MC-Data" in Table El- 
These differences were estimated by propagating the dif- 
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FIG. 6: (color online). One-dimensional scans of — 2 A In/! as a function of magnitudes (left) and phases (right) of the 
resonances /o(980), /o(1710), /2(2010), and XcO (top to bottom). The horizontal dashed lines mark the one- and two-standard 
deviation levels. 



ferences, in the control sample, between background- 
subtracted data and signal MC, into the fit PDFs. 

For the branching fraction measurement, we assign a 
systematic uncertainty due to the error on the calcula- 



tion of N BB ("AW in Table [V} and to the K s recon- 
struction efficiency. We correct the Kg reconstruction 
efficiency by the difference between the efficiency found 
in a dedicated K s data sample and that found in simu- 
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TABLE III: The interference fractions FF(k,j) among the intermediate decay amplitudes for Solution 1. Note that the diagonal 
elements are those defined in Eq. (|28l) and detailed in Table ITT1 The lower diagonal elements are omitted since the matrix is 
symmetric. 





fo(9SO)K u s 


/o(1710)Jfg 


/ 2 (2010)J^ 


NR 




/o(980)J3 


0.44 


0.07 


-0.02 


-0.80 


0.01 


/o(1710)Jfg 




0.07 


-0.01 


-0.17 


-0.0003 


/ 2 (2010)^§ 






0.09 


0.02 


0.0002 


NR 








2.16 


-0.02 


XcoK° s 










0.07 



TABLE IV: Summary of measurements of branching fractions 
(£?). The quoted numbers are obtained by multiplying the 
corresponding fit fraction from Solution 1 by the measured 
inclusive B° — > KgKgKg branching fraction. The first un- 
certainty is statistical, the second is systematic, and the third 
represents the signal DP-model dependence. 



Mode 




B 


[xlO" 6 ] 


Inclusive B° - 


+ K%K%K% 


6.19 ±0.4 


8 ±0.15 ±0.12 



/o(980)JC§, / (980) ^ tfgtfg 
fo(17W)K° s , /o(1710) 
/a(2010)X§, / a (2010) 
NR, K%K%K% 
XcoK° s , xco ^K° S K° S 



■ +1-3 



K U S K U S 0.50 ZqZ ± 0-04 ± 0.10 
K%K° S 0.54 toil ± 0.03 ± 0.52 

+2.2 



0.46 ±R'?5 ± 0.02 ± 0.21 



lation. We assign the uncertainty on the correction as a 
systematic error ("Kg reco" in Table lv| . 



categories of the tag-flavor assignment. Furthermore the 
decay rate is convolved with the per-event At resolution 
7?. s ig(At, a At), which is described by the sum of three 
Gaussians and depends on At and its error a At- For an 
event i with tag flavor qt ag , one has 



£ , si g (A£,a-At;g'tag,c) 

-\At\/r E 

1 + qta 



(33) 



4r B o 
+<Zta g (-D) f 



AD C 



S sin( Arrid At) — C cos(AnidAt) 



® U^^t, a At), 

where qt ag is defined to be ±1 (—1) for B tag = B (B tag = 
B ), tqo is the mean B° lifetime, and Amj is the mixing 
frequency 16]. The widths of the B° and the B° are 
assumed to be the same. 



B. Event selection and backgrounds 



IV. TIME-DEPENDENT ANALYSIS 

In Sec. IIV Al we describe the proper-time distribution 
used to extract the time-dependent CP asymmetries. In 
Sec. IIV Bl we explain the selection requirements used to 
obtain the signal candidates and suppress backgrounds. 
In Sec. IIV Cl we describe the fit method and the approach 
used to account for experimental effects. In Sec. IIV Dl we 
present the results of the fit and finally, in Sec. IIV El we 
discuss systematic uncertainties in the results. 



A. Proper-time distribution 

The time-dependent CP asymmetries are functions of 
the proper-time difference At — tcp — t tag between a 
fully reconstructed B° — >• K S K S K S decay (B C p) and 
the other B meson decay in the event (-Btag), which is 
partially reconstructed. The observed decay rate is the 
physical decay rate modified to include tagging imperfec- 
tions, namely (D) c and AD C ; the former is the rate of 
correctly assigning the flavor of the B meson, averaged 
over B° and B°, and the latter is the difference between 
D c for B° and B°. The index c denotes different quality 



We reconstruct B° — > KgKgKg candidates either 
from three Kg — > 7r + 7r~ candidates, or from two Kg — > 
7r + 7r~ and one K s —> tt°tt , where the tt° candidates 
are formed from pairs of photons. The vertex fit require- 
ments are the same as in the amplitude analysis, and also 
the requirement that the charged pions of at least one of 
the K s have hits in the two inner layers of the vertex 
tracker. The K° s candidates in the 5° 3K s (ir + 7r~) 
submode must have mass within 12MeV/c 2 of the nom- 
inal K° mass [l6| and decay length with respect to the 
B vertex between 0.2 and 40 cm. In addition, combi- 
natorial background is suppressed in both submodes by 
imposing that the angle between the momentum vector 
of each K s (tt + tt^) candidate and the vector connect- 
ing the beamspot and the K s (tt + tt~) vertex is smaller 
than 0.2 radians. Each K s decaying to charged pions 
in the B° —> 2K s (ir + ir~)K s (Tr Tr a ) submode is required 
to have decay length between 0.15 and 60 cm and 7r + 7r~ 
invariant mass less than 1 1 MeV from the world aver- 
age Kg mass pj. The K% decaying to neutral pions 
in the B° ->• 2A^(7r+7r")i^g(7r 7r ) submode must have 
7T 7T° invariant mass between 0.48 and 0.52GeV/c 2 . Ad- 
ditionally, the neutral pions are selected if they have 77 
invariant mass between 0.100 and 0.141 GeV/c 2 and if 
the photons have energies greater than 50 MeV in the 
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TABLE V: Summary of systematic uncertainties. The model uncertainty is dominated by the variation of the line shapes due 
to the contribution of the poorly measured /2(2010). 



Parameter Fit bias B-bkg Discr. Vars MC-Data N B g Kg reco Sum Model 



B(B° -> K° S K° S K Q S ) [HT 6 ] 


0.011 


0.030 


0.053 


0.015 


0.067 0.111 0.145 


0.120 


FF /o(980) 


0.013 


0.056 


0.006 


0.001 


0.058 


0.190 


FF /o(1710) 


0.007 


0.001 


0.001 


0.001 


0.007 


0.016 


FF / 2 (2010) 


0.005 


0.001 


0.003 


0.001 


0.006 


0.084 


FF NR 


0.024 


0.083 


0.023 


0.001 


0.090 


0.344 


FF X c0 


0.002 


0.000 


0.001 


0.000 


0.002 


0.034 


Ph [rad] /o(980) 


0.008 


0.018 


0.014 


0.000 


0.024 


0.177 


Ph [rad] /o(1710) 


0.011 


0.020 


0.001 


0.003 


0.023 


0.185 


Ph [rad] / 2 (2010) 


0.044 


0.014 


0.004 


0.002 


0.046 


0.684 


Ph [rad] Xco 


0.039 


0.011 


0.010 


0.007 


0.042 


0.498 



laboratory frame and a lateral energy deposition pro- 
file in the electromagnetic calorimeter consistent with 
that expected for an electromagnetic shower (lateral mo- 
ment [24| less than 0.55). The fact that we do not model 
any PDF using sideband data allows a loose requirement 
on Toes and AE in the time-dependent analysis, namely, 
5.22 < toes < 5.29GeV/c 2 and -0.18 < AE < 0.12 GeV. 
In case of multiple candidates passing the selection, we 
proceed in the same way as in the amplitude analysis. We 
use the same NN as in the amplitude analysis to suppress 
continuum background. 

With the above selection criteria, we obtain signal 
reconstruction efficiencies of 6.7% and 3.1% for the 
B° -> 3K° s (7r+n-) and B° -> 2K° s (Tr+7r')K g(vr 7r ) sub- 
modes, respectively. These efficiencies are determined 
from a DP-model MC sample generated using the results 
of the amplitude analysis. We estimate from MC that 
2.1% of the selected signal events are misreconstructed 
for B° -> 3ifg(7T + 7r"), while the figure is 2.4% in 
B" -> 2ifg(7r + 7r-)^(7r°7r°), and we do not treat these 
events differently from correctly reconstructed events. 
Because of the looser requirements, there are more back- 
ground events from B decays than in the amplitude 
analysis, in particular in the B a -> 2Kg(-K + it~)Kg('K a Ti a ) 
submode. These backgrounds are included in the fit 
model and are summarized in Table IVI1 As the analysis 
is phase-space integrated, we cannot model the Xco res- 
onance separately, and its contribution to the CP asym- 
metries could cloud deviations in the charmless contri- 
butions. We therefore apply a veto around the invariant 
mass of this charmonium state. 



C. The maximum-likelihood fit 

We perform an unbinned extended maximum- 
likelihood fit to extract the B n -> KgK^K^ event yields 
along with the S and C parameters of the time-dependent 
analysis. 

The fit uses as variables rn-ES) AE, the NN output, At, 
and <7At- The selected on-resonance data sample is as- 
sumed to consist of signal, continuum background, and 
backgrounds from B decays. Wrongly reconstructed sig- 



nal events are not considered separately. The likelihood 
function d for event i is the sum 

£ * =^N j V l j (m ES ,AE,At,(T± t ,NN-,q ta z,c,p), (34) 

3 

where j stands for the species (signal, continuum back- 
ground, one for each B background category) and Nj is 
the corresponding yield; (/tag, c, and p are the tag flavor, 
the tagging category, and the physics category, respec- 
tively. 

To determine <7tae and c we use the B flavor-tagging 
algorithm of Ref. [25J. This algorithm combines several 
different signatures, such as charges, momenta, and de- 
cay angles of charged particles in the event to achieve 
optimal separation between the two B flavors. This pro- 
duces six mutually exclusive tagging categories. We also 
retain untagged events in a seventh category; although 
these events do not contribute to the measurement of 
the time-dependent CP asymmetry they do provide ad- 
ditional sensitivity for the measurement of direct CP vi- 
olation [26j. 

The two physics categories correspond to 
B° -> 3A^(7r+7r-) and B° -> 2K° s (TT+7r-)K° s (tt°tt°) 
decays. The PDF for species j evaluated for event % is 
given by the product of individual PDFs: 

fj( m ES: AE, At, (TAt, NN; 9tagj c, p) = (35) 
Vj(m ES ;p) V](AE;p) 7>j(NN;c,p) V]{At,a At ;q tag ,c,p). 

To take into account the different reconstruction of the 
two submodes, we use separate PDFs for the two physics 
categories. Separate NN and At PDFs are included for 
each tagging category within each physics category. The 
separate At PDFs for the two physics categories allow us 
to fit the S and C parameters either separately for the 
two submodes, or together. The total likelihood is given 
by 

C = exj>(-Y f N j )'[[U (36) 

3 i 
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TABLE VI: Summary of B background modes included in the fit model of the time-dependent analysis. The expected number 
of events takes into account the branching fractions (B) and efficiencies. In case there is no measurement, the branching fraction 
of an isospin-related channel is used. All the fixed yields are varied by ±100% for systematic uncertainties. 



Submode 






Background mode 


Varied 


B [xl0" b ] 


Number of events 


JD — r Ol\ g [TV TV 


1 




l<~0 lt-0 JfO 


no 


O A 


U. / 1 








J\ 5-A- <ji\ 


no 


z/.o 


y.oo 








K° K° K+ 


no 


1 1 5 


A 97 






B° - 


-V -fnpntvfll cptipHp Hppavsl 


yes 


nnt annlipflblp 


21.7 






B+- 


+ {charged generic decays} 


yes 


not applicable 


15.5 


B" -> 2K%{-k+-k- 






K%K%Kl 


no 


2.4 


0.67 








K° S K° S K* 


no 


27.5 


5.3 








K° S K° L K* 


no 


27.5 


0.3 








K° S K° S K+ 


no 


11.5 


2.9 








K%K%K*+ 


no 


27.5 


7.2 






B°- 


-> {neutral generic decays} 


yes 


not applicable 


73.6 






B+- 


■> {charged generic decays} 


yes 


not applicable 


73.8 



1. At PDFs 

The signal PDF for At is given in Eq. (|33|) . Param- 
eters that depend solely on the tag side of the events 
(namely, (D) c and AD C ) are taken from the analysis of 
B —J- ccK^ decays [27|. On the other hand, parame- 
ters that depend on the signal-side reconstruction, due 
to the absence of direct tracks from the B decay, can- 
not be taken from modes that include such direct tracks. 
This is the case for the parameters that describe the res- 
olution function, which are found in a fit to simulated 
events. A systematic uncertainty for data-MC differences 
is assigned using the control sample B° — > J/tfjKg, as ex- 
plained in Sec. IIVEI 

For continuum events we use a zero-lifetime compo- 
nent. This parametrization is convolved with the same 
resolution function as for signal, with different parame- 
ters that are varied in the fit to data. The parameters 
of this PDF are not separated in the tagging categories. 
The small contribution from e + e~ — > cc events is well 
described by the tails of the resolution function. For B 
background events we use the signal PDF, with resolution 
parameters from the BABAR B — > ccK^ analysis [27l |. 
The parameters S and C are set to zero and varied to 
assign a systematic uncertainty. 



2. Description of the other variables 

The toes aud AE distributions of signal events are 
parametrized by an asymmetric Gaussian with power- 
law tails, as given in Eq. (|24|) . and, for toes, a small 
additional component, parameterized by an ARGUS 
shape function [2(J, to correctly describe misrecon- 
structed events. The means in these two PDFs for 
B° —> 3Kg(ir + ir~) events are allowed to vary in the fit 
to data, and the other parameters are taken from MC 
simulation. For the NN distributions of signal we use his- 



togram PDFs taken from MC simulation for each physics 
and tagging category. The Toes, AE, and NN PDFs for 
continuum events are parametrized by an ARGUS shape 
function, a straight line and the sum of power functions 
from Eq. (|25p . respectively. All continuum parameters, 
except for C2 and C3 of the NN PDF, are allowed to vary 
in the fit. All the fixed parameters are varied, within the 
uncertainties found in a fit to sideband data, to estimate 
systematic errors. 

All the B background PDFs are described by fixed 
histograms taken from MC simulation. 



D. Results 



The 



maximum-likelihood fit of 3261 candidates in 
ir~) submode and 7209 candidates in 
the B° -> 2K^(TT + TT-)K^(n Tr ) submode results in the 
event yields detailed in Table IVTlI 



the B° ->■ 3if § 



TABLE VII: Event yields for the different event species, 
resulting from the maximum-likelihood fit for the time- 
dependent analysis. U B + B~ (B°B°) bkg" represents back- 
ground from charged (neutral) B decays. Quoted uncertain- 
ties are statistical only. 



Species 


3K° s (ir+n-) 




Signal 
Continuum 
B+B- bkg 
B°B° bkg 


201 tit 
3086 t 5 5 l 

01 _ 24 
n+31 
a -30 


62 til 
7086 til 

45+34 
40 -30 
4 +38 
4 -29 



The fit result for the time-dependent CP-violation pa- 
rameters S and C is 



s = -0.94 ±8;|* 

C = -0.17 ±0.18, 
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where the uncertainties are statistical only. The correla- 
tion between S and C is —0.16. We use the fit result to 
create /Plots of the signal distributions of At, the time- 
dependent asymmetry, and the discriminating variables. 
Figure [7] shows the At s Vlots for the combined fit result 
and for the individual submodes. Figure |8] shows the sig- 
nal distributions and Fig. |H] the continuum background 
distributions of the discriminating variables. The distri- 
butions shown in these three figures illustrate the good 
agreement between the data and the fit model. 

We scan the statistical-only likelihood of the S param- 
eter for both submodes and for the combined fit. The 
result, on the left-hand side of Fig. [101 shows a sizable 
difference between the S values for the two submodes; 
the level of consistency, conservatively estimated from 
the sum of the two individual likelihood scans, is approxi- 
mately 2.6er (a p- value of 1.0% with 2 degrees of freedom). 
This value is obtained including only the dominant statis- 
tical uncertainty and neglecting the small correlation be- 
tween the CP-violation parameters. The results obtained 
when S and C are allowed to vary individually for each 
of the submodes are S = -1.42 ±g;|J, C = -0.14 for 
B° -> 3K° s {it+it-) and S = 0.40 C = 0.19 t° H for 
B° -> 2K s (ir + ir-)K s (TT a Tr a ). In both cases the quoted 
uncertainties are statistical only. 

As there is some correlation between the S and C pa- 
rameters, we perform a two-dimensional statistical like- 
lihood scan of the combined likelihood, which is then 
convolved by the systematic uncertainties on S and C 
(systematic uncertainties are discussed in Sec. lIVEl ) The 
result is shown on the right-hand side of Fig. [TU] We find 
that CP conservation is excluded at 3.8 standard devia- 
tions, and thus, for the first time, we measure an evidence 
of CP violation in B° -> K S K%K% decays. The differ- 
ence between our result and that from P° — > ccK^** is 
less than 2 standard deviations. The scan also shows that 
the result is close to the physical boundary, given by the 
constraint S 2 + C 2 < 1. 



E. Systematic uncertainties 

The systematic uncertainties are summarized in Ta- 
ble I Villi The "MC s tat" uncertainty accounts for the 
limited size of the simulated data samples used to cre- 
ate the PDFs. The "P rcC o" uncertainty propagates the 
experimental uncertainty in the measurement of tag- 
side-related quantities taken from (27j to our measure- 
ment. The "P-bkg" contribution results from the un- 
certainty in the CP content and the branching fractions 
of fixed yields in the model of background from B de- 
cays. The dominant "MC-Data: At" systematic uncer- 
tainty is due to possible differences between data and 
simulation concerning the procedure used to obtain the 
signal B decay vertex from tracks originating from K§ 
decays. We quantify this uncertainty using the control 
sample P° -> J/ipK° s by comparing the difference be- 
tween At values obtained with and without the J/0 in 



data and simulation. We then propagate the observed 
differences and their uncertainties to the resolution func- 
tion. We use this new resolution function to refit the 
data and obtain an estimate of the effect on S and C. 
We also use the samples B° — > J/ipK s(7T + 7r _ ) and B° — > 
J/ip Kg(TT°Tr ) to estimate simulation-data differences for 
the other variables in the submodes B° — > 3Kg(ir + Tr~) 
and B° -s- 2X°(7r+7r-)^(7r 7r°), respectively. This con- 
tribution is referred to as "MC-Data: Discr. Vars" in Ta- 
ble I Villi The "Fit Bias" uncertainty is evaluated using 
fits to fully reconstructed MC samples. It accounts for 
effects from wrongly reconstructed events and correla- 
tions between fit variables. The "Vetoes" uncertainty is 
related to the veto on the invariant mass. It is evaluated 
using events that pass the veto in pseudo-experiments 
studies. Finally, the "Misc" uncertainty includes contri- 
butions from doubly-Cabibbo-suppressed decays, silicon 
vertex tracker alignment, and the uncertainties in the 
boost of the T(45). These contributions are taken from 
the BABAR B -> caFr w analysis [27|. 



V. SUMMARY 

We have performed the first amplitude analysis of 
B° —> KgKgKg events and measured the total inclusive 
branching fraction to be (6. 19 ±0.48 ±0.15 ±0.12) x 10~ 6 , 
where the first uncertainty is statistical, the second 
is systematic, and the third represents the signal DP- 
model dependence. We have identified the dominant 
contributions to the DP to be from / (980), / (1710), 
/2(2010), and a nonresonant component, and measured 
the individual fit fractions and phases of each compo- 
nent. We do not observe any significant contribution 
from the so-called /x(1500) resonance seen in, for exam- 
ple, B + -t K + K~K + 0. The peak in the invariant 
mass between 1.5 and 1.6GeV/c 2 can be described by 
the interference between the /o(1710) resonance and the 
nonresonant component. We see some hints from the 
/ 2 (1525) and /o(1500) resonances that could also con- 
tribute to this structure, but due to limited sample size 
we cannot make a significant statement. Future inves- 
tigations of the KK system could shed more light on 
the situation. Furthermore we have performed an up- 
date of the phase-space-integrated time-dependent anal- 
ysis of the same decay mode, using B° — > 3Kg(TT + n~) 
and B Q -> 2i^(7T+7r-)i<:g(7r 7r ) decays, with the final 
BABAR dataset. We measure the CP-violation parameters 
to be S = -0.94 i°;^± 0.06 andC = -0.17± 0.18± 0.04, 
where the first quoted uncertainty is statistical and the 
second is systematic. These measured values are con- 
sistent with and supersede those reported in Ref. Q. 
They are compatible within two standard deviations 
with those measured in tree-dominated modes such as 
B a — » J/ipKg, as expected in the SM. For the first time, 
we report evidence of CP violation in B° — > KgK^Kg 
decays; CP conservation is excluded at 3.8 standard de- 
viations including systematic uncertainties. 
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FIG. 7: (color online). Signal sVlots (points with error bars) and PDFs (histograms) of At (left) and the derived asymmetry 
(right) for the B° -¥ 3A'g(7r+7r") submode (top), the B° -> 2K^(tt + 7v~)K^(7v°tt°) submode (middle), and the combined fit 
(bottom) . In the At distributions on the left-hand side, points marked with x and solid lines correspond to decays where Stag 
is a B° meson; points marked with o and dashed lines correspond to decays where B tag is a B° meson. Points of the asymmetry 
s Vlots that are outside the range of a figure are marked by arrows. 



TABLE VIII: Summary of systematic uncertainties on the S and C parameters. 



Source 



MC s tat 0.002 0.001 

Breco 0.004 0.003 

B-bkg 0.032 0.012 

MC-Data: At 0.045 0.027 

MC-Data: Discr. Vars 0.021 0.004 

Fit Bias 0.022 0.018 

Vetoes 0.006 0.004 

Misc 0.004 0.015 

"Sum 0.064 0.038 
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FIG. 8: (color online). Signal a Vlots (points with error bars) and PDFs (histograms) of the discriminating variables: mBS (top), 
AE (middle), and the NN output (bottom) for the B° -> 3K° s {tt+-k-) submode (left) and for the B° -> 2ifg(7r + 7r-)A"g(7r 7r°) 
submode (right). Below each bin are shown the residuals, normalized in error units. The horizontal dotted and full lines mark 
the one- and two-standard deviation levels, respectively. 
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FIG. 9: (color online). Continuum a Vlats (points with error bars) and PDFs (histograms) of 7ties, AE, the NN output, and 
At (top to bottom). Plots on the left-hand side correspond to the B° — > 3Kg(-K + n~) submode, and on the right-hand side to 
the B° -> 2K° s (-k + ty- )Kg(iY iv ) submode. In the At distributions, points marked with x and solid lines correspond to decays 
where -Stag is a B° meson; points marked with o and dashed lines correspond to decays where -Btag is a B° meson. Below each 
bin are shown the residuals, normalized in error units. The horizontal dotted and full lines mark the one- and two-standard 
deviation levels, respectively. 
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FIG. 10: (color online). One-dimensional statistical scan of — 2Aln£ as a function of S (left) and the two-dimensional scan, 
including systematic uncertainty, as a function of S and C (right). In the left-hand plot, red points marked with x correspond 
to the B° — > 3Kg(-ir + -K~) submode, blue points marked with o to the B° — > 2K%(-k + -k~)K%(-r -k°) submode, and black points 
marked with * to the combined fit. In the right-hand plot, the gray scale is given in units of V~2A ln£. The result of the 
BABAR analyses of B a — > ccK^ decays [53] is indicated as a white ellipse and the physical boundary (S 2 + C 2 < 1) is marked 
as a gray line. The scan appears to be trimmed on the lower left since the PDF becomes negative outside the physical region 
(i.e., the white region does not indicate that the scan flattens out at 5a). 
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